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ABSTRACT 

The partial differential eguation which describes steady-state diffusion 
from an infinite line source has been replaced with a set of simultaneous ordi- 
nary differential equations solved on an electronic analog computer. One space 
dimension^ distance downwind^ was represented by computer time; the other, 
height, was replaced with finite differences,. Solutions were obtained for 
constant, power law, and logarithmic wind profiles, and for diffusion of partic- 
ulates which can settle out and deposit on the ground. 

All solutions were obtained with one basic computing circuit „ Each prob- 
lem required only a particular setting of the coefficient potentiometers in the 
circuit. Implementation of this circuit required only 9 integrating amplifiers 
and 26 coefficient potentiometers, available in any medium sized computer. 

The solution's accuracy was measured by comparing the computer plots with 
the analytical solution for constant wind profile. This measured the total 
error due to the finite difference approximation and to computer errors. The 
solution's accuracy was found to be 5^ or better over most of the field. 
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1. INTRODUCTION 

Diffusion due to atmospheric turbulence can be described by the parabolic 
diffusion equation,, This is the direct physical approach in contrast to that 
which invokes statistical concepts. Through the use of such an equation and 
its boundary conditions, diffusion problems are conveniently specified, but 
analytical solutions are not so readily obtained. So great is the difficulty, 
indeed, that it is unusual to find a solution that fits problems of interest. 
With the aid of various computing devices, such solutions can be found „ Given 
such a device, in this case an electronic analog computer, it is of interest 
to explore the form of the solutions as a function of the various equation 
parameters. 

The physical problem considered here is that of steady-state diffusion 
from an infinite line source, oriented normal to the wind„ The solution vlll 
specify the concentration of the diffusing substance as a function of position 
In the region after the source has been emitting at a steady rate for a long 
time„ It -will be obtained as a plot of concentration as a function of dis- 
tance downwind for discrete height Intervals, As examples of this process, 
consider diffusion of exhaust gases from automobiles on a busy road, or a uni- 
form, slowly advancing grass fire. The physical model, that is, the idealized 
case. Is shown in Fig. 1 and consists of a line source of infinite extent along 
the y axis. The x axis is oriented downwind which is nonmal to the source and 
the z axis extends vertically. The wind vector is assumed to have no components 
in the y or z directions. Since the source has no height or width, the concen- 
tration at the source must be infinite, 
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2. THE PARABOLIC DIFFUSION EQUATION 



The general diffusion equation is 



dt 



1 



^^ Sx 



Sy 



K ^ 






K. 



^ 

Sz 



(1) 



where 



dt ^t Sx hj dz 



and 



X = concentration of the diffusing substance 

t = time 
Kj_ = coefficient of eddy diffusivity in the i direction 

u = mean wind speed in the x direction 

V = mean wind speed in the y direction 

w = mean wind speed in the z direction. 
In accordance with the physical model^, the derivatives with respect to y are 
zero, and the wind components v and w are zero. The steady-state provision 
eliminates the partial derivative with respect to time. If we also neglect 
the diffusion term in x compared to the translation term, that is, if 



u 



Bx 



» 



Sx 



-^ Sx, 



then the diffusion equation becomes 



T7 In 



u 



Sx h'z. 



dz. 



(2) 



With the equation in this form, the solution will depend upon the fomi of u 



and Kg and upon the boundary conditions. 
If u and K^ are given by 



m 



u(z) = Ui( — j , KgCz) 



Ki 



n 



and the boundary and continuity conditions are 



KgSx/Sz -^0 as z->0, x>0 



^ ^ 00 along X = z = 



/ ux(x, z) dz = Q for all x > 
o 



(3) 



where Q is the rate of emission of material per unit length of source, then the 
solution given by Sutton is 



x(x, z) 



uir(s) 



u-i 



(m - n + 2) Kix 



exp 



Uiz"^-^+^ 



(m - n + 2) Kix 



(4) 



where 



s = 



m + 1 
m - n + 2 



x(x,z) is the mean concentration of the diffusing substance. This is a very 
useful solution since it can be used for any power law wind profile provided 
that m - n + 2 > 0. For the very important case where the wind profile is 
given by u oc log z, K^ °^ z, no analytical solution has been found. 

An approximate solution for the logarithmic wind profile, u « log z, 

2 
K «: z, was given by Karplus and Allder. They used the electric analog which 



is similar in principle to the technique reported here. It consisted of building 



a ladder network of resistors and capacitors which satisfied a set of simul- 
taneous differential equations which was an approximation to the partial dif- 
ferential equation (2) . 

The electronic analog approach is different primarily in that a standard 
electronic analog computer is utilized instead of biiilding a special circuit 
as was done by Karplus and Allder. This has the advantage of using standard, 
general purpose equipment, and is much more flexible. 



3. VARIOUS MATHEMATICAL-PHYSICAL MODELS 

The basic mathematical model which expresses the physical system shown 
in Fig. 1 is Eq, (2) , The special cases to be considered here are 

(a) Constant velocity profile with standard boundary conditions [those 
specified in (j) ]; u = const^ K^ = const, 

(b) Constant velocity profile with standard boundary conditions and 
gravitational settling; u = const ^ K^ = const, w = fall speed. 

(c) Constant velocity profile with ground absorption; u = const, 

^0, 



Kg, = const, -^ 

^" z=0 
(d) Constant velocity profile with gravitational settling and ground 



absorption; u = const, Kg = const, w = fall speed, ^ 



Sz 



^0. 
z=0 

(e) Power law wind profile with standard boundary conditions; 

u(z) = ui (^) , Kz(z) = Ki (^) 

(f) Logarithmic wind profile with standard boundary conditions; 
u(z) = ux log(-^) , KgCz) = ^i(~)' 



The complete model Is 



^1 Sx h 
A( z) OX dz 



B(z) '^] ^Z^ 



w^ (5) 

dz 



Sz, 

where the functions A(z) and B(z) represent the wind speed and eddy duffusivity 
terms, respectively. Thus 

u(z) = -^ , K^Cz) = KiB(z) 
A(z) 

Settling of particulate matter is represented as a component of the wind directed 



earthward. The iDoundary conditions for (5) are 



alona; x = z = 



/ ux(x, z)dz 
o 



for all X > 






as z -> 0, X > 



for cases a^ 'b, e, and f. 



Y ^ ^ at z = 0, X > 

^ Bz 

for cases c and d. 



To facilitate handling Eq, (5) and to increase the utility of the solu- 
tions, it will be expressed in nondimensional terms. Set S = x/Xn ^^'^ ^ ~ ^/^o 
where the reference parameters x,o ^'^^ ^o '^iH ^^ defined later. Then Eq.. (5) 
may be written 



I = '^'^ I 



B(Z) 






+ W :r— 
SZ 



(6) 



where 



X = 



W = 



Kn 



X 



ZqUi 



ZqW 



Ki 



and the symbol X is reserved for incorporation of some more constants which 
will appear later. 



k. PREPAHATION FOR ANALOG COMPUTi\TION 

Since this problem is "being prepared for solution on an analog computer^ 
we must observe the restriction that the analog computer can integrate only 
with respect to one independent variable, namely, computer time. As stated 
[Eg. (6) ], the model has two independent variables x and Z so that one of 
them, in this case Z, must be eliminated. This is done by replacing deriva- 
tives with respect to Z by finite differences as follows: 



SS 
SZ 



3n+i ■" '=n-i 



az 



SS 
SZ 



n 



1 


1 + (AZ) 1 

2 '■'■'2 

Sn+1 - Sn 
(AZ) . 


Sn " Sn_i 


(AZ)^ 


(AZ) 1 

n-2 - 



Implementation of finite differences requires a new model to clarify the meaning 
of the position n and the increment AZ. Such a model is presented in Fig. 2 
which shows the atmosphere divided into 10 layers in the vertical. The thick- 
ness of the layers increases exponentially with height so that the thickness of 
the top layer is infinite. The reason for this exponential distribution is 
discussed more fully later. Each layer has a thickness AZ, is infinite in the 
y direction, and semi- infinite in the x direction. All properties of each 
layer are lumped in the center of the layer so that the model may be thought 
of as a stack of planes spaced exponentially. Each plane incorporates the 
properties of the layer of atmosphere which it replaces. For convenience, the 
planes are labeled 1 through 10 and represent the comrputlng stations. Thus the 
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Fig. 2. Finite difference model showing the vertical coordinate 
axis relations. 
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notation Sj^(x) represents the concentration in the nth layer as a function 

of the distance downwind. 

The first partial derivative with respect to Z required for the settling 

term must he centered at the station n and must span from n - 1 to n + 1 so 

that the increment of Z is Z^+i - Z^-i. If the layer spacing were uniform, 

this would he just 2(AZ)o The notation ( AZ) ^ indicates the span from Zv, to 

n+^ 

Z„ . T , and ( AZ) „ indicates Z -r - Z t. 

Wow the meaning of Xq ^^'^ ^n '^^^ ^® defined in terms of the finite dif- 
ference model. In (3) , the continuity condition 

CO 

Q = / u x(x,z)dz (x > 0) 
o 

relates the source strength Q to the concentration and wind speed „ In the 
finite difference model, this expression holds at x IS 0, so at x = 0, 

a _ 
Q = / u X dz 
o 

= a u X 

since the concentration is taken to be constant throughout the lowest layer, 
which extends to a height a = 0„1053Zq, and is zero above this layer. Labeling 
the concentration at x = as Xn^ then 



V = _:2 (7) 

^o 0.1053ZqU ^" 

If u is a function of height, the value used here is the average over the height 
interval g z g 0„1053Zq. 

In the model the source is a region of uniform emission through some height 

11 



increment and the spacing of the layers is set accordingly as a function of 
this height increment. 

Using finite difference expressions for the derivatives with respect to Z, 
Eq. (6) will be replaced with a set of 10 simultaneous differential equations, 
one for each layer. The form of these equations is 



dS 
dX 



An 



n 



(AZ) 



n 



^^ (AZ)n+i " n--| (AZ) 1 

2 2 



^ (AZ)n+i + (AZ)^_i ^ 



The new independent variable X incorporates any arbitrary constants needed. 
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5. MINIMIZATION OF ERROR DUE TO FINITE DIFFERENCE APPROXIMATION 

Since Eq. (8) is only an approximation to Eq. (6) , it is important to 
analyze the inherent error and to find vays of minimizing it. Considering 
the case where the layer spacing in the vertical is uniform^ the finite dif- 
ference approximation for the second derivative of S with respect to Z is 
given by 



a^s 


1 


Sn+i - Sn 


Sn - Sn-i 


Bz2 


. " AZ 


AZ 


AZ 



Sn+i " 2Sn + Sn-i 



To evaluate the error involved^ expand Sj^+i and S^-i about the point n using 
a Taylor expansion,^ 



Sn+i - S; 



'n 



AZ SS 
i: BZ 



+ (AZ)"5^S 
21 SZ2 



^ (AZ)-' Sfs 
3i ' 5Z3 



n 



(9) 



Sn-i 



S - ^ ^ 

" " i: Sz 



+ Uz)f afs 

2.' SZ^ 



( AZ) ^ afs 
31 SZ^ 



n 



(10) 



Add the expressions (9) and (lO) and solve for B^S/SZ^, 



SZ2 



n 



Sn+i - 2Sn + Sn-i 
(AZ)2 



+ e 



where 



12 aZ4 



n 



(11) 



The error that we are seeking is e and it is a function of the station spacing 
and of the gradient at a given station. The error can be reduced by increasing 
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the number of stations vMch decreases the height increment. However, since 
the amount of equipment involved is proportional to the number of stations, 
we must keep the number to a minimum. Since the gradient is greatest near 
the source and falls off with distance from the source, the contribution of 
the gradient term to the error will decrease with height. Thus we space the 
stations closer together near the ground with small height increments and let 
the height increments increase exponentially with height. This reduces the 
error for a given number of stations and provides more data in the lower layers 
where the greatest changes occur. 
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6. EQUATIONS FOR VAEIOUS MATHEMATICAL-PHYSICAL MODELS 

Equation (8) with the appropriate hotindary conditions for each of the 
seven cases is given below. 

(a) Constant velocity profile. The simultaneous differential equations are 



dS 
dX 



100( AZ) 



=n+i 



^n ^n " Sn-i 



(AZ) T 
n+^ 



(AZ) 



1 
n-- 



(12) 



where 



X 


= 


100 X 




Si 


= 


1.00 


at X = 


So 


= 


Si 


for all X 


Sio 


:= 





for all X 



An arbitrary factor, in this case 100, is introduced in each case to facilitate 
handling the equations on the computer. Also, since the independent variable X 
corresponds to computer time, the arbitrary factor controls the solution speed 
which must be within the limits set by the frequency response of the two-coordinate 
plotter, used to record the solution, on the one hand, and by the patience of 
the operator on the other. 

The boundary conditions must be restated in terms of the finite differences. 
The source strength is no longer Infinite at a line but finite at an area as 
given by (7) . Thus, the concentration ratio Si in the first layer is unity 
at X = 0. The requirement that the flux across the boundary, represented by 
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the ground at n = — ^ be zero is satisfied by setting the concentration at 
Sq, a virtual station below ground, equal to that at Si for all values of X. 
The flux at n = — is given by 



as 

az 



^^Z-^ = 



AZ 

2 



The concentration in the tenth layer must be zero because that layer is in- 
finite in extent. Since S^q = for all values of X, there are only 9 active 
cells and 9 simultaneous differential equations when ( 12) is expanded. 
(b) Constant velocity profile with gravitational settling. 



dS 
dX 



n 



Sn+i - Sn Sn - S; 



n-i 



lOO(AZ)^ L (AZ)^^i (AZ)_i 



n+i ' ' n-^ 



1 Sn+i - Sn-i 



100 (AZ) . + (AZ) 1 
— n+— n— = 

2 2 2 2 



(13) 



X = 100 X 
Si = 1.00 at X = 
So = Si for all X 

Sio = for all X 

W = 1 for all X and Z. 

This equation is different from the previous one only in that the settling term 
has been added. Since the purpose here is only to show how to incorporate set- 
tling and the nature of the effect on the solutions, the settling coefficient 
was arbitrarily set at W = 1. If the eddy diffuslvlty Ki were 10^ cm^ sec""^ 
and Zq were 10 cm, the fall speed would be 10 cm sec"'^ since W = wZq/Ki. 
(c) Constant velocity profile with ground absorption. 
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dS 
dX 



1 



n 



3n+i 



- s. 



'n ^n 



Sn - S- 



'n-1 



lOO(AZ) L (AZ) . 



^^ 



(AZ)^_, 

2 



(1^) 



vhere 



X = 100 X 
Si = 1.00 
So = 7S1 



5lO 



= 



at X = 

for all X (7 = 0.5) 

for all X. 



Equation (ik) is just like (l2) with the exception of one of the "boundary condi- 
tions. If we define an absorption coefficient 7 such that when 7=0 there is 
complete absorption, the ground acts as a sink. When 7 = 1.00^ there is complete 
reflection at the ground. This model cannot distinguish between passage through 
the boundary and deposition on it. For this case, 7 = 0.5^ which means that 
the concentration in the virtual layer n = is Just one-half the concentration 
in the first layer for all X. In the case of deposition, enough material is 
deposited on the boundary to provide the required concentration in the virtual 
layer. 

(d) Constant velocity profile with gr avlt at i onal settling and ground 
absoi-ption. 



dS 
dX 



100( AZ) 
n n 



(AZ) ^1 (AZ) 1 
n+— n--- • 
2 2 



1 ^n+x ~ Sn-i 

100 (AZ) 1 + (AZ) 1 
n+-| n-| 



(15) 



where 



X = 100 X 



Si = 1.00 



at X = 
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So = 7S1 for all X (7 = O.5) 
Si5 = for all X. 

This case provides the logical combination of the two previous ones. The ma- 
terial is assiiraed to have a finite fall speed and is allowed to deposit on the 
boundary surface. 

(e) Power law wind profile. 



6 6 



dS 
dX 



where 



^(~) \(-) ^Sn+i - Sn) - (^) (Sn - Sn-i)l (16) 



n+-^ n— „ 

2 2 



X = 86.36 -^ ^ X 
Si = 1.00 at X = 

So = Si for all X 

Sio = for all X. 

In comparing the distribution of concentration in various wind profile regimes^ 
differences occur due to a change in the average wind speed and to the wind 
speed distribution. If we set Zq = Zj_ to eliminate the arbitrary constants 
raised to the ^ power in X, comparison of this case to, say case (a) or (f), 
at equal values of X will not show both effects but only that due to a dif- 
ferent wind distribution. This is so because the number 86.56/IOO compensates 
for the different average values of K and u. Comparison at equal values of x 
will show both effects. 

(f) Logarithmic wind profile. 



dS 
dX 



63.72 
n 



(log Z) 
(AZ) 



-1 1 



n 



(I) 1 ^'-'^ - '-^ - (I) 1 ^^- - ^--^] 

n+i n-| 

(IT) 
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where 

X = 63.72 X (Zi = Zo) 

Si = 1.00 at X = 

So - Si for all X 



Sio - for all X 



log Z = natural logarithm of Z. 
As in the previous case, compensation has been introduced for the average 

u and K. 
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7. MALOG MODEL 

The physical model now consists of a stack of planes spaced exponentially 
In the z direction. It is supposed that each plane has a capacity for contain- 
ing some of the diffusing suhstance just equal to that of the layer of atmosphere 
that it replaces. At a given point in x,, a traverse in the z direction shows 
the properties of the atmosphere lumped at discrete intervals. This model 
can he simulated with an electric analog., where each of the circuit components 
is analogous to some property of the atmosphere. A segment of such a circuit 
is shown in Fig. J. The analogies shown were developed from the similarity 
of the node differential equations for a simple case^ 



dV 

Cn — 

dt 



(Az)n u ^ 
dx 



= ^ (Vn+i - 2Vn + Vq-i) 
n ^ 



(Xn+i - 2xn + Xn-i) 



n \"^/n 



(Az). 

To perform, the simulation^ set the resistor Rn proportional to Az/k, the capaci- 
tor Cn proportional to Az u and voltage V proportional to concentration. The 
problem time t will represent distance downstream. .All these relations are 
stated as proportionalities so that scale factors can be used to provide reason- 
able values of resistance and capacitance. In general, one may relate electri- 
cal resistance to atmospheric resistance to diffusion and electrical capacitance 
to the capacity of a layer of the atmosphere to hold the diffusing substance. 
The complete network analog is shown in Fig. h along with the electronic analog 
and a representation of the physical model. 
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Fig, k. Electronic analog, passive network analog, and physical, 
finite difference model of the atmosphere. Note that each group 
of components represents a layer of the atmosphere. 
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The passive circuit analog was not used in tMs study, but it is helpful 
to construct it on paper because of the insight it provides. This is particu- 
larly useful in detemiining the method of applying the boundary conditions. 
To insure that Si(0) = 1.00, the Si capacitor is disconnected from the circuit 
and charged to an arbitrary voltage which is taken to be unityo At the position 
(computer time) X = 0, the capacitor is switched into the circuit. The re- 
quirement that Sq = Si for all X is simply provided by making the resistor 
Rq which connects the virtual station Sq to Si infinite. Then the flux (cur- 
rent) between them is zero as required. If we want Sq = (nonref lee ting 
boundary) , connect resistor Rg to ground as shown in Fig. k. However, the 
requirement that Sq = 0,5 Si is difficult to implement in the passive analog. 
Again, to insure that Siq = 0, make the position 10 at ground potential. 

The electronic analog circuit could be designed directly from the pas- 
sive network or from the set of simultaneous differential equations. The lat- 
ter method is much the easier course after Eqs. (l2)-(l7) are reduced to the 
form 



dS 
dX 



= ^n^n+i - "bn^n + ^n^n-i - (l8) 

n 



The principal element in the electronic analog is the integrating amplifier 
which can be used to represent one station as shown in Fig. 5- I^ "the passive 
network analog, there was a direct correspondence between the circuit elements 
and the physical problem, while in the electronic analog, there exists a one- 
to-one correspondence between the circuit elements and the defining equation (l8) . 
The coefficients a, b, and c are set on the correspondingly designated coefficient 
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HIGH-GAIN, d-c 
AMPLIFIER 



o- 



COEFFICIENT POTENTIOMETER 



Fig. 5- Typical station in the active or electronic analog circuit. 
The concentration is represented as a voltage^ and equation coef- 
ficients as potentiometer settings,, 
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potentiometers. The integrating amplifier sums the three inputs, integrates 
the sum, and inverts the sign. Ordinarily, the gain of the integrating ampli- 
fier is unity, set by the 1 megohm input resistor and 1 microfarad feedback 
capacitoro If one of the coefficients is greater than unity, the correspond- 
ing input resistor used is 0,1 megohm which provides a gain of 10, This is 
necessary because the maximim setting of the coefficient potentiometers is 1. 

The electronic circuit is much more flexible than the passive network 
permitting almost any conceivable form of boundary conditions to be set up 
readily. 
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8. COMPARISON OF MALOG MD MALYTICAL SOLUTIONS 

The form of the solution is a voltage at the output of each amplifier which 
is proportional to the concentration at the appropriate level in Z and which 
is a continuous function of X. When recorded on a two-coordinate plotter^ the 
result is a family of curves representing the concentration at discrete levels 
of Z. 

The evaluation of errors in this tyj)e of computation is especially im- 
portant since there are two very different possible sources. One is the approxi- 
mation involved in using finite differences and the other is due to the com- 
puter. Previously, we indicated how the error due to the use of finite differ- 
ences could he minimized hut said nothing ahout the magnitude of the errors. 
The procedure here will he to compare the computer solution of case (a) con- 
stant wind profile, with the analytical solution which may he obtained from 
Eq, {k) by setting m = n = 0. 



uir(|) 



Ui 



1 

2 



exp 



2 



UiZ 



(19) 



.kKiX 

Equation (19) is the solution for diffusion from a line source of infinite 
concentration while the computer solution, is for a vertical area source of 
limited height and of finite concentration. There seems to be some problem 
In comparing these solutions. The actual source strengths can be related by 
substituting Eq. (7) into Eq. (19) and, for convenience, setting 
Ki = ui = Zq = 1. 
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S = -^ = °^i251_exp ' "' 



z 



(20) 



In terms of the nondimensional parameters X and Z, (20) is 



S = 0.59^1X exp 



2 _._. I 25Z' 



2' 



X 



(20a) 



¥e can determine whether the finite difference solution is a good approxima- 
tion to (20a) by comparing the form of the solutions with respect to X and Z 
and by comparing values of S at selected points in the X,Z space. A simple 
relation between S and X exists at Z = such that 



= O.59J+IX ^ . (21) 

Z=0 



Then we can observe the relation of S to the exponential term by finding the 

,2 



slope of logarithm of S with respect to Z 



Slope = log Si - log Sg ^ _ 25 (22) 

Z? - Zf X 



Table I shows how well (2l) and (22) are satisfied. The exponential term seems 
to fit well out to X > 100 where the slope was very small. At X & 200^ all 
the concentrations are small, S g .OU so that a constant error which is small 
for S becomes a large percentage error when S is small. This is shown in 
Table II where computer values of concentration are compared with analytical 
values . 

The most significant feature of Table II is that the magnitude of the 
error is fairly constant. The percentage error increases as the concentration 
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TABLE I 

DETEMINATIOW OF THE FIT OF THE FINITE DIFFERENCE SOLUTION 
TO THE ANALYTICAL SOLUTION 



X 


1 
^^^ Z == 


-X [slope] 


1 


0.5^9 


2U.8 


5 


0.586 


26.3 


10 


0.573 


25.6 


25 


0.572 


25.3 


50 


0.566 


25.5 


100 


0.573 


25.4 


200 


Oo572 


32.1 


300 


0.532 


43 = 2 


Analytical 






Values 


0.5941 


25.0 



TABLE II 

CX)MPAEISON..DF ANALYTICAL AND: FINITE DIFFERENCE SOLUTIONS. 

AT SELECTED POINTS 



n 



Concentration 
Analytical 



Magnitude 
of Error 



Percentage 
Error 



1 
2 

5 

k 

5 



1 
2 

5 

h 

5 
6 

7 



25 



100 



0.119 


0.0050 


4.2 


Oell6 


0, oo4o 


5.5 


0.109 


0.0040 


3.7 


0.0987 


0. 0045 


4.6 


0„083i 


0,0038 


4.6 


0.0628 


0.0027 


h,3 


0„039^ 


0.0019 


4.8 


0.0174 


0.0011 


6.3 


0.00325 


0.00025 


7.7 


0.0594 


0.0021 


3.5 


0.0590 


0.0020 


3A 


0.0582 


0,0021 


3,6 


0.0567 


0.0020 


3.5 


0.0543 


0.0020 


3 = 7 


0.0502 


0.0015 


3.0 


0.0447 


0.0013 


2.9 


0.0368 


0.0016 


h,3 


0.0242 


0,0022 


9=1 
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decreases Indicating that there is a limit to useful computation. When the 
concentration ratio S falls to less than 0.001, computer error in the form 
of noise becomes significant. Virtually all the error shown in Table II is 
due to the limitations of the finite difference model. The error in this 
computation is less than 5^ provided that S > 0,001 and except in layer 9° 
Greater errors should be expected in this layer since it is the largest 
active one and is next to the infinite layer for which zero concentration is 
assumed. 
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9. ADDITIONAL RESULTS OF COMPUTATION 

To present the information more compactly^ the solutions have heen re- 
plotted on log-linear paper as shown in Fig. 6. This presents concentration 
as a continuous function of height although what was actually obtained was a 
histogram as shown in Fig. 7- 

Taking Fig. 6 as a standard for comparison^ it is seen from Fig. 8 that 
the effect of ground absorption is to produce lower concentrations in the pro- 
files. Figure 8 also shows that the maxiLmum concentration found in the profile 
for X > occurs above the ground, thus indicating that the plume axis lifts 
off the ground. Figure 9 shows that gravitational settling produces profiles 
of somewhat lower concentrations than in Fig. 6 and with more slope. In this 
case, of course, the plijme does not rise. The combination of settling and 
absorption or deposition is shown in Fig. 10. 

Solutions due to power law wind profile and logarithmic profile are shown 
in Figs. 11 and 12, The level of concentration in the profiles is very nearly 
the same but the slopes are slightly different. In both cases, the profiles 
for small values of X approximate to straight lines. 

To give an intuitive appreciation for the nature of the solutions obtained, 
they are presented in Figs. I3-I8 as lines of constant concentration in the 
X, Z space. 
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Fig. ?• Histogram of concentration versus height. Constant wind 
profile at X = 50. This figure indicates the real form of the 
computer solutions. 
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10. SUWMAEY MD CONCLUSIONS 

Solutions to the diffusion equation for the case of the infinite line 
source have "been found for three types of wind speed profile and for the case 
of diffusion of particulates with an appreciable fall speed. These are sample 
problems that were chosen to demonstrate the method. In the problems selected, 
the functions of wind speed and eddy diffusivity with height were specified 
in analytical form and then translated to discrete steps in height. Since 
they are used in this form, one could as easily start with functions expressed 
graphically or In any arbitrary manner. Indeed, using the passive network 
analogy as a guide, it is possible to design the computing circuit directly 
from the physical problem. Even the layer spacing used need not follow some 
mathematical rule; it could be chosen to fit the problem. 

The wind speed and eddy diffusivity were assumed constant with distance 
downwind. Changing them with respect to x would involve changing resistors 
in time which is not difficult to implement with servo multipliers which are 
usually found in analog computer installations. In the same way, the boundary 
conditions could be made a function of x. Thus one could simulate diffusion 
over undulating ground. 

The accuracy obtained here was better than 5fo over most of the field. 
This could be improved to 1^ by increasing the number of cells inasmuch as 
analog computers are available with component quality of O.Olfo of full scale. 
Extension beyond ifo solution accuracy is not practical with this technique, 
nor is there any theoretical or practical need for g2:-eater accuracy than this. 

^5 



The advantages of the analog technique in problems of this type are the 
speed and convenience with which a solution may he obtained and the ease with 
which atmospheric parameters may be varied. In addition to its value as a re- 
search tool^ the analog method should be a valuable teaching aid in the analysis 
of atmospheric diffusion. Taking advantage of the direct relation between 
mathematical operations and computer components, an instructor could demonstrate 
before a class, the effects of changing parameters on the solution obtained. 

This method may be extended to problems of finite line sources, area 
sources both horizontal and vertical, and to elevated point sources. 
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APPENDIX 

The following tables list the equation coefficients in the form of Eq. 
(18) . These are the coefficients that were applied to the computer in the 
form of coefficient potentiometer settings. Equation (18) is repeated for 
convenience: 



dS 
dX 



- S'n^n+i ~ ^nSn "^ ^nSn-i 



n 



(18) 



TABLE A- I 



CONSTANT WEMD PROFILE 



n 


a 


b 


c 


1 


0.8571 


1 

0.8571 





2 


0.6780 


l.kkk 


0.7661 


3 


0.5226 


1.089 


0.5665 


k 


0.3887 


0.8il22 


0.^51+5 


5 


0.2735 


0.6011 


0.3276 


6 


0.1782 


0.i|015 


0.2233 


7 


O.lOji^ 


0. 21+16 


0.1382 


8 


O.OI+833 


0.1218 


0.073^8 


9 





O.Oij-131 


0.02820 



TABLE A- II 



CONSTANT WIND PROFILE WITH GROUND ABSORPTION 



Same as Table A-I except bx = I.308 for 7 = 0.5 
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TABLE A- III 



CONSTANT WIND PROFILE WITH GRAVITATIONAL SETTLING 



n 


a 


b 


c 


1 


0.9046 


0.8571 





2 


0.7205 


l.U4i^ 


0.7236 


3 


Oo5595 


1.089 


0.5292 


k 


0,ij-202 


0.81|-22 


0,4220 


5 


0.5009 


0.6011 


0.3002 


6 


0.2006 


0.ii015 


0.2007 


7 


0.1208 


0.2416 


0,1208 


8 


0.06068 


0.1218 


0,06113 


9 





0.04131 


0,02099 



TABLE A- IV 

CONSTANT WIND PROFILE WITH GRAVITATIONAL SETTLING 
AND GROUND ABSORPTION 



Same as Table A- III except bi = I.308 for 7 = O.5 



TABLE A-V 



POWER LAW WIND PROFILE 



n 


a 


b 


c 


1 


0.1956 


0.1956 





2 


0,2499 


0,3982 


0.1483 


3 


0.2663 


0.4690 


0.2027 


4 


0.2534 


0.4721 


0.2187 


5 


0.2213 


0.4255 


0.2042 


6 


0.1759 


0.3494 


0.1735 


7 


0.1240 


0,2551 


0,1311 


8 


0,07137 


0.1549 


0.08355 


9 





0.06504 


0.03984 
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TABLE A- VI 



LOGARITHMIC WI1\ID PROFILE 



n 


a 


b 


c 


1 


O.IU19 


0.l4l9 





2 


0.2220 


0.3403 


0.1183 


5 


0.2550 


0,4369 


0.1819 


h 


0.2523 


0.4589 


0.2066 


5 


0.2246 


0.4230 


0,1984 


6 


0.1795 


0.3497 


0.1702 


T 


0.1264 


0.2548- 


0.1284 


8 


0.06782 


0.1450 


0.07719 


9 





0.06286 


0.03774 
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